# ------------------------------------------------------------------------------
# Create 31-365d mace/death outcome regression tbls for tested v untested
# Author: Cassidy Shubatt <cshubatt@gmail.com>
# To run: bash 04_regress_longterm_adverse.sh
# ------------------------------------------------------------------------------

# Libraries --------------------------------------------------------------------
library(here)
library(yaml)
library(data.table)
library(tidyverse)
library(glue)
library(lfe)
library(stargazer) # latex tables

u <- modules::use(here("lib", "util.R"))
temp <- here(
  "code", "03_score_diagnostics", "temp",
  "04_regress_longterm_adverse"
)
if (!dir.exists(temp)) {
  dir.create(temp)
}
overnight_lab <- ""

# Load Data --------------------------------------------------------------------
message("Loading data...")
paths <- read_yaml(here("lib", "filepaths.yml"))
longterm_outcomes <- readRDS(paths$analysis$longterm_outcomes)
cohort <- readRDS(glue(paths$analysis$test_cohort)) %>%
  filter(!exclude) %>%
  u$safe_left_join(longterm_outcomes) %>%
  mutate(
    yhat = p__ensemble__stent_or_cabg_010_day__tested,
    yhat_quintile = tile_stent_or_cabg_005_tested,
    yhat_quartile = tile_stent_or_cabg_004_tested
  )

# Outcome Regressions ----------------------------------------------------------
message("Getting regression table for longterm outcomes...")
adverse_outcomes <- c(
  "macetrop_or_death_31_to_365", "macetrop_31_to_365_pos", "death_365_day"
)
t_vars <- c("test_010_day")
risk_vars <- c("yhat", "yhat_quartile")

reg_df <- cohort
for (outcome in adverse_outcomes) {
  reg_df[[outcome]] <- reg_df[[outcome]] * 100
}

for (t_var in t_vars) {
  for (risk in risk_vars) {
    reg_list <- list()
    reg_list_nointer <- list()

    # Regressions
    for (outcome in adverse_outcomes) {
      message(
        "Regressing ", outcome, " ~ ", risk, " + ", t_var, " + ",
        risk, ":", t_var, " in "
      )
      form <- glue("{outcome} ~ {risk}*{t_var} | 0 | 0 | ptid") %>% formula()
      fit <- felm(form, data = reg_df)
      reg_list <- c(reg_list, list(fit))

      form_nointer <- glue("{outcome} ~ {risk} + {t_var} | 0 | 0 | ptid") %>% formula()
      fit_nointer <- lm(form_nointer, data = reg_df)
      reg_list_nointer <- c(reg_list_nointer, list(fit_nointer))
    }

    # Stargazer Tables
    omits <- c("f", "ser", "ll", "aic")
    reg_tbl <- stargazer(
      reg_list,
      omit.stat = omits, no.space = TRUE, dep.var.caption = ""
    )
    save_fp <- file.path(temp, glue("longterm_regs__adverse__{risk}__{t_var}.tex"))
    write(reg_tbl, save_fp, append = FALSE)

    reg_tbl_nointer <- stargazer(
      reg_list_nointer,
      omit.stat = omits, no.space = TRUE, dep.var.caption = ""
    )
    save_fp <- file.path(temp, glue("longterm_regs__adverse_interacted__{risk}__{t_var}.tex"))
    write(reg_tbl_nointer, save_fp, append = FALSE)
  }
}

message("Done.")
